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INTRODUCTION 


This  presentation  discusses  the  development  of  a  technique  to  numerically  model  the  reaction  of 
composite  materials  that  are  undergoing  microdamage.  Damage  on  the  microstructure  of  composites, 
such  as  microcracks  and  fracture,  can  have  a  significanteffect  the  global  response  of  the  system.  The  most 
notable  effect  is  in  the  change  of  material  propertiesofthe  composite  as  microcracks  form  and  propagate. 
Though  empirical  data  can  be  used  to  approximate  the  material  properties  for  an  ideal  non-damaged 
specimen  [1],  performingexperi  mentation  on  damaged  compositesand  extracting  a  correlation  between 
microcrack  size  and  or  quantity  to  material  properties  can  be  extremely  difficult. 

One  approach  to  extract  material  properties  from  a  composite  material  without  having  to  perform 
experimental  testing  is  a  numerical  approximation  method,  called  the  homogenization  method  [2].  This 
method  uses  asymptotic  expansions  of  field  variables  about  macroscopic  values  and  provides  overall 
effective  properties  as  well  as  microscopic  stress  and  strain  values.  A  limitation  with  this  method, 
however,  is  that  it  assumes  uniformity  of  the  macroscopic  fields  within  each  representative  volume 
element  (RVE).  Hence,  this  method  breaks  down  in  critical  regions  of  high  gradients  such  as  cracks.  This 
further  complicates  its  use  for  modeling  microcrack  phenomena,  as  it  is  in  the  regionsof  global  cracks  that 
microcracking  will  be  formed. 

To  overcome  the  drawbacks  of  the  existing  methods  [3,4]  proposed  a  structural  based  enrichment 
method  based  on  the  principles  of  partition  of  unity  which  allows  enrichment  of  the  approximation  space 
in  localized  sub  domains  usingspecialized  functions  that  may  be  generated  based  on  a  priori  information 
regardi  ng  asymptotic  expansionsof  local  stress  fields  and  microstructure.  One  of  the  major  advantages  of 
the  partition  of  unity-based  enrichment  strategies  is  that  the  enriched  local  function  space  maybe  easily 
varied  from  one  node  to  the  other. 

In  this  presentation  we  review  the  development  of  RVE  models  which  will  correlate  composite  material 
properties,  based  on  microcracking,  and  their  incorporation  into  a  modified  structural  based  enrichment 
method  to  capture  the  microdamage  effects.  This  approach  can  be  used  in  regions  where  macrocracks 
are  assumed  to  have  initiated.  In  the  next  section,  we  review  the  homogenization  and  enrichment 
methods  and  discuss  the  implementation  of  damaged  microstructures  into  the  enrichment  technique.  In 
section  3,  we  present  an  example  which  demonstrates  the  effectiveness  of  our  approach. 

EXPERIMENTATION 

We  have  divided  this  section  into  three  parts.  Section  2.1  and  2.2  review  the  concepts  of  the 
homogenization  method  and  structural  enrichment  technique,  respectively.  Section  2.3  discusses  the 
implementation  of  damaged  microstructures  into  the  structural  enrichment  technique. 
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Figure  1.  Multiscale  Homogenization  Approach 
Review  of  the  Homogenization  Approach 

The  homogenization  theory  [2],  assumes  that  a  function,  composed  of  multiple  scales,  can  be  designated 
as 

fY(x)  =  f{x,y(xj)  ^ 

wherexandy  represent  a  macroscale  and  microscale,  respectfully  (see  Figure  l).The  scales  can  be  related 
by  a  scale  factor,  y,  such  that 

x 

y  =  y  U1 


A  key  feature  of  the  homogenization  approach  is  the  assumption  of  local-periodicity  on  the  microscale. 
Thus,  if  we  have  a  representative  volume  element  (RVE),  as  shown  in  Figure  1,  with  a  side  length  of  Y we 
can  write 

f(xi,yJ)  =  f(xi,yJ+kYj)  [3] 

where  /cis  the  periodicinterval.  The  partial  derivative forthe  function  with  respecttoxcan  be  simplified 
to 

d*t/(Wy)  =  dxJ  +  Y^dy.f  Mi 


Using  asymptotic  expansion,  we  can  write  the  displacement,  stain  and  stress  fields,  with  respect  to  the 
scale  factory,  as: 
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u\(x,y)  =  u°i(x,y)  + yu}(x,y)+ 0(y2)  [5a] 

ejj  O,  y)  =  y-1  £[/  (x,  y)  +  e?  (x,  y)  +  ye}  (x,  y)  +  0  (y2 )  [5b] 

<*ij(x.y)  =  Y~1(Jij1(x>y)  +  °i°(x,y)  +yai1(x,y')  +  0(y2)  [5c] 

where 

ffij  ~  C ijkl£kl  [5] 


is  the  constitutive  relationship  for  the  constituents  of  the  RVE.  In  this  presentation  we  will  assume  that 
the  constitutive  behavior  of  each  constituent  of  the  RVE  is  linear  elastic  and  known. 


To  determine  the  effects  of  the  microstructure  across  the  domain,  we  begi  n  by  substituti  ng  the  asymptotic 
expansions  into  the  equilibrium  equation  and  rearranging  the  terms  to  give  us: 


'2  \{Cim£kl)  y.  +  y 

+y°  +  (Cijkl£kl)  —  fi 

L  i  >yj 


1  ( Cijkl£ki L)  +  i^ijkl£kl) 

L  ’xj  .yj 


+  0(y)  =  0 


[7] 


In  equation  [7],  isthe  elastic  coefficients  and  isthe  body  force.  By  settingeach  ordertermto  zero 
we  derive  the  following  equations 

u°(*,y)  =  u°U)  [8] 

[C,mOUm„  +  ^sM')]yi=  0  e£2RVE  PI 

[CmAn,]  ~ftB  =  0  &  [10] 

,Xj 

Equation  [8]  states  that  is  only  a  function  of  the  macroscale.  Equation  [9]  is  the  governing  equation  for 
the  microscale  RVE,  used  to  derive  the  effective  elastic  coefficients,  where  ip^nki 's  the  symmetric  part 
of  the  spatial  gradient  of  H^nk  defined  as 


mnkl  2 


LjS 

mnl,yk  '  r2mnk,yi > 


[ii] 


hffy)  is  a  Y-periodic  tensor  that  relates  the  macroscale  to  the  first  order  displacement  perturbation  term 
u}(x,y) ,  such  that: 


uj(x,y)  =  H^ni(y)e^nx(x)  [12] 

Equation  [10]  is  the  governing  equation  for  the  macroscale.  In  the  equation,  C  represents  the  effective 
elastic  coefficients  and  fB  is  the  effective  body  force.  They  are  respectively  defined  as: 
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Cijkl  —  yRVE  f{}RVE Cijl<lOklmn  +  %rinkl)  3flRVE 

?iB  =  - fjB  anRVE 
VRVE  Jn RVE  1 


[13] 

[14] 


Structural  Enrichment  Method 

Within  this  section  we  review  the  partition  of  unity  paradigm  and  how  it  is  i  mplemented,  followed  by  the 
derivation  of  the  enrichment  functions. 


Interpolation  Functions 

The  partition  of  unity  paradigm  requires  that  the  sum  of  all  functions  at  a  given  point  in  a  domain  is  equal 
to  one, such  that 

N 

Yj(p?M  =  1  [15] 

i= 1 


In  the  above  equation,  <p;°  is  the  Ith  partition  of  unity  function  and  N  is  the  total  number  of  partition  of 
unity  functions  within  the  domain,  O.  Thus,  the  sum  of  all  the  functions  at  a  given  point,  x,  within  O  will 
equal  1. 

Based  on  this  concept  we  define  the  global  approximation  space,  for  a  method  using  interpolation 
functions,  as: 

Vh,p  =  YN  (p?V,h’p  c  i/Hn)  [16] 

*—‘i= i 


where  the  superscripts  h  and  p  are  the  size  of  the  element  and  the  polynomial  order,  respectively.  H1  is 
the  first  order  Hilbert  space  and  is  the  local  approximation  space  at  /,  defined  as: 

V^p  =  spanm^{pm(x)}  c  H^B, DH)  [17] 

In  equation  [17],  ^is  an  index  set,  B/  is  the  support  for  node  /  and  pm(x)  is  composed  of  polynomials  or 
other  functions.  From  the  above  equations  we  can  write  any  function  vhp  within  the  global  approximation 
space  as: 

N 

VKp(x)  =  II  hlm(x)aIm,  vh’p  6  Vh’p  [18] 

/  =1  met 


where 


[191 

is  the  interpolation  function  at  node  /,  corresponding  to  the  mth  degree  of  freedom,  and  aIm  is  the 
associated  degree  of  freedom. 
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For  FEA,  the  standard  shape  function,  Nh  which  can  be  referenced  in  numerous  sources,  satisfies  the 
partition  of  unity  criteria  for  cpj.  For  a  typical  finite  element  analysis,  correlating  the  definition  of  a  finite 
element  shape  function  to  the  interpolation  function  in  equation  [19],  requires  that  pm  —  span{  1}. 
Reference  [5]  discusses  other  partition  of  unity  functions  and  basis's  that  can  be  implemented. 


Enrichment  Functions 


In  the  previous  section  we  defined  a  function  in  the  global  approximation  space  as  the  sum  of  the  partition 
of  unity  functions  multiplied  by  a  group  of  basis  functions  and  associated  degrees  of  freedom.  Focusing 
on  the  basis  functions,  p„[xj,  we  state  that  if  an  arbitrary  function  is  included  in  the  local  basis  of  all  nodes 
within  the  domain  and  assuming  aIm  —  Smn  V/  then 

N 

^  0/°  O) 

1  =  1 


2^  2^  0/°  =  Pn  (*)  I 


/=  1  mcf 


Equation  [20]  demonstrates  that  it  is  possible  to  exactly  reproduce  a  function  overthe  entire  domain.  This 
fundamental  conceptenablesthe  consistency  of  the  interpolation  functionsto  be  adjusted  based  on  the 
population  of  the  basis  functions.  Thus,  we  can  expend  the  basis  function  to  include  relevant  functionsto 
enforce  a  known  displacementfield.  Recallingthatthe  basisfora  standard  FEA  analysis  ispm  =  span{  1}, 
and  incorporating  the  asymptotic  response  about  a  crack  tip,  a  common  enrichment  technique  is  to  adjust 
the  basis  to 

pm(x)  =  span{  l.yfr)  [21] 

Using  this  adjusted  basis  within  the  region  of  a  crack  tip  will  enable  the  interpolation  functions  to 
reproduce  the  Vr  field.  If  the  crack  tip  is  notat  the  origin  of  r  the  degree  of  freedom  associated  with  the 
enrichment  function  should  solve  to  zero,  thus  removing  the  effects  of  the  enrichment  term  on  the 
displacementfield. 

To  derive  the  structural  based  enrichmentfunction,  we  begin  byexaminingthe  asymptotic  expansion  of 
the  displacementfield  given  in  equation  [5a].  Substituting  equations  [8]  and  [12]  into  [5a]  and  assuming 
0(y2)  approximation  gives  us: 

(x, y )  =  u°i  (x)  +  y  (y)s^nx(x))  [22] 

Fora  single-scale  problem,  we  define  the  displacementfield,  assuming  linear  consistency  in  2D,  as 

N 

u\  (x)  =  u°  O)  =  ^  ^  0°  (x)pm (x)  aImi  [23] 

1  =  1  me^ 


where  <p° and  pm  are  /V,  and  pm  =  span{  1}. 

To  perform  a  multiscale  analysis  we  need  to  incorporate  the  second  term  on  the  right  hand  side  of 
equation  [22],  We  focus  on  the  Y-periodic  tensor,  Hs(y).  Recal  ling  that  y  is  a  function  of  x,  equation  [2]  and 
assuming  hffyfx))  is  the  dominant  term,  we  can  write  the  relation: 
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N 

Y  (Hmni  (y)£%mx  (*))  =  ^  <Pi  CK*))/Wii 

7=1 


[24] 


Substituting  [23]  and  [24]  into  [22]  gives  us: 

N  N 

u\  (x)  =  II  0/°O)PmO)  ^ Imi  ♦ 1  0/°(*)^mni(y(*))^  I  mni  [25] 

1=1  me^  1=1 

Rearranging  [25]  and  incorporating  the  degrees  of  freedom  vector  6  into  a  reduces  [25]  to 

N 

u\  O)  =  ^  ^  0°  U)Pmi  (*)  almi  t20] 

7=1 

The  enrichment  function  can  be  defined  as: 

Tinl(x)  =  Hgnl(y(x))  [27] 

One  significant  observation  is  that  incorporating  the  enrichment  function  requires  that  the  basis  be 
dependent  on  the  coordinate  direction.  For  example,  a  two  dimensional  finite  element  analysis  would 
have  the  enriched  basis  in  the  Xj-direction,  as 

Pmi  =  span{l,  tffn  (y(x)),//2s21  (yW),  tfmCyU))}  t28] 

while  for  the  x2-direction  the  enriched  basis  is 

Pm2  =  span{l,  H(12  {y{xj),H$22  (yW)»  Hi 122  (yW)}  [29] 

In  genericform  fora  finite  element  simulation,  the  basis  can  be  written  as 

pmi  =  span{  1,  T^ni  (x)}  [30] 

Implementation  the  Multiscale  Enrichment  Technique 

The  approach  to  implementing  structural  based  enrichment  varies  depending  on  the  governing  method. 
In  this  presentation  we  will  focus  in  the  FEA  approach.  Reference  [3]  gives  complete  details  on  the 
implementation  of  the  multiscale  enrichment  into  FEA.  In  this  section  we  will  briefly  review  how  the 
algorithm  is  applied. 

The  analysis  is  based  on  a  homogenization-like  integration  (HLI)  scheme  [3].  The  sol'n  requires  the 
enrichment  function,  calculated  using  equation  [9],  and  the  material  properties,  which  are  only  known  on 
the  microscale,  y.  To  implement  this  technique  and  resolve  the  aforementioned  issue,  integration  will  be 
performed  on  an  RVE  transformed  to  the  macroscale,  x,  coordinate  system.  The  HLI  scheme,  shown  in 
Figure  2,  is  based  on  the  assumption  that  there  is  a  significant  deference  between  the  two  scalesand  that 
the  Gauss  points  used  for  FEA  are  sparse  with  respect  to  the  size  of  the  RVE.  In  this  approach,  an  RVE  is 
centered  on  each  Gauss  point  within  macroscale  problem.  An  integrand,  /,  on  the  macroscale  is  replaced 
with  the  following 
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nguass  nguass 

1  =  I  Jvdn=  ^  wijiv(xfp)=  £  wjiV^j_^Hx)dnf™  [31] 

a  i  =  1  i  =  l  ni 

Where  ir  is  an  arbitrary  function,  xfp  is  the  /"'Gauss  point,  nfVE  is  the  domain  of  an  RVE  centered  on  the 
Fh  Gauss  point,  J  is  the  Jacobian  and  l/l/,isa  weight  associated  with  the  ?h  Gauss  point.  [3]  has  demonstrated 
that  the  accuracy  of  HU  improves  with  decreasing  size  of  the  RVE.  Thus,  this  implementation  works  well 
when  the  difference  in  scales  is  significant.  However,  as  the  scales  become  closer  and  the  therefore  the 
volume  of  the  RVE  increases,  equation  [31]  will  develop  convergence  errors. 


Figure  2.  Implementation  of  multiscale  enrichment  into  FEA 

Flow  Chart  of  Structural  Based  Enrichment  Algorithm 

In  this  section,  we  discuss  the  algorithm  for  the  structural  based  enrichment  method  depicted  in  Figure  3. 
We  begin  the  analysis  by  deriving  the  Hsfrom  the  RVE  and  material  properties  provided.^  is  then  inputted 
into  the  macroscale  structural  problem  which  solves  forthe  displacement,  stress  and  strain  fields. 


Figure  3.  Flow  Chart  of  Structural  Enrichment  Approach 
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Implementation  of  Damaged  Micros  tincture 

We  now  focus,  in  this  section,  on  the  RVE and  how  microdamage  can  be  incorporated  into  the  model.  For 
many  applications  the  material  used  in  the  multiscale  model  is  some  type  of  fiber  composite,  such  that 
the  microstructure  can  be  presented  as  a  fiberwithin  a  matrix  (see  Figure  4a)  fora  2D  application,  and  a 
section  of  the  fiber  weave  (see  Figure  4b)  for  3D  simulations. 


Figure  4.  Microstructures  for  a  Fiber  Composite 

For  normal  conditions,  these  RVEs  can  be  used  to  represent  the  microstructure;  however,  these  multiscale 
methods  are  specifically  designed  for  critical  regions,  such  as  in  the  vicinity  of  a  global  crack,  where  fiber- 
matrix  is  likely  to  be  damaged  on  the  microscale  as  well. 


Figure  5.  Damaged  Matrix  Microstructures  for  a  Fiber  Composite 

To  account  for  microdamage,  the  RVEs  are  modified  to  incorporate  varying  levelsof  cracks.  Figure  5 gives 
examples  of  several  cases  where  cracks  are  imbedded  into  the  microstructure  matrix. 

The  cracks  on  the  RVE  are  modeled  as  geometric  damage  in  the  simulation,  such  that  there  is  no 
connectivity  across  a  crack  on  elements  within  the  RVE  neigh  boring  the  crack.  Asa  pre-step  to  solving  the 
global  problem,  we  solve  equation  [9]  on  the  damaged  RVE  to  produce  a  set  of  damaged  enrichment 
functions,  H^-.  The  damaged  enrichmentfunctionsare  then  incorporated  into  model  withinthe  vicinity 
of  a  global  crack,  see  Figure  6. 
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Figure  6.  Critical  Region  Around  Global  Crack 

The  algorithm  for  the  modified  enrichment  method  is  given  pictorially  in  Figure  7.  We  begin  the  analysis 
by  deriving  the  H5  and  Ff'0  from  the  RVEand  material  properties  provided.  f/5andi/5_Dare  then  inputted 
into  the  macroscale  structural  problem,  in  the  form  of  enrichment  functions.  The  enrichment  functions 
derived  from  the  damaged  RVE,  HS~D,  are  appliedto  elementswithin  the  vicinity  of  the  critical  region, 
while  the  remaining  elements  within  the  global  model  are  enriched  with  functions  derived  from  the 
standard  RVE,  Hs.  The  profile  of  the  critical  region  is  dependent  on  criteria  provided  by  the  modeler. 


Figure  7.  Critical  Region  Around  Global  Crack 


RESULTS 

We  examine  a  two-dimensional  cross  section  of  a  damaged  composite  wafer,  depicted  in  Figure  8, 
containing  a  center  crack.  A  pressure  load  is  applied  to  the  upper  and  lower  surfaces  to  force  open  the 
crack  and  a  temperature  change  is  implemented  to  induce  thermal  stress.  The  additional  enrichment 
functions  used  to  capture  the  thermal  stress  is  reviewed  in  [5].  In  this  example  we  assume  plane  stress 
conditions  and  examine  the  upper  right  quadrant  of  the  cross  sectional  area.  We  assume  the  composite 
is  composed  of  graphite  fibers,  all  orientated  perpendicularto  the  cross  section,  within  an  epoxy  matrix. 
The  material  properties  are  given  in  Table  1.  The  RVEs,  are  modeled  using  500  quadratic  2nd  order 
elements,  where  the  enrichment  functions  generated  from  the  damaged  RVE's  are  applied  to  the  global 
region  encompassed  by  the  yellow  boundary  and  the  enrichment  function  gene  rated  from  the  standard 
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RVE  isappliedthroughoutthe  remainderof  the  model. The  boundary  conditions  on  the  global  model  are 
a  40  MPa  pressure  pul  ling  on  the  upper  and  lower  surfaces  of  the  wafer,  and  a  heat  source  increasing  the 
upper  and  lower  surfaces  to  1080°C.  The  results  are  compared  to  a  fine  mesh  composed  of  1.8e6 
quadratic  2nd  order  elements  modeled  down  to  the  microscale. 


Table  1.  Material  Properties  of  RVE 


Material 

Elastic  Modulus 
(MPa) 

Poisson's 

Ratio 

Thermal  Conductivity 
(mW/(mmC)) 

Coefficient  of  Thermal 
Expansion  (C1) 

Epoxy 

3450 

0.35 

4.47 

5.4xl0'5 

Graphite 

700 

0.485 

2.4 

1.2  xl0‘5 

40MPa 


28mm ' 
“20mm 


1028C 


0.06mm 


Figure  8.  Damaged  Composite  Sample 


Figure  9.  RVE  Damaged  Orientations 

We  begin  the  analysis  by  examining  several  crack  orientations  shown  in  Figure  9.  For  this  analysis,  the 
homogenization  simulation  is  run  using  900  quadratic  2nd  order  elements,  with  the  effective  material 
properties  generated  from  a  damaged  RVE  applied  to  the  critical  region  around  the  crack.  To  compare  the 
orientations,  we  examine  the  strain  energy  resulting  from  the  simulation  with  respect  to  the  results  from 
the  fine  mesh  which  produced  a  strain  energy  of  882.4  MPa.  Table  2  shows  the  resulting  relative  error  in 
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the  strain  energies.  The  results  show  that  the  orientation  does  not  make  significant  difference  in  error 
generated.  The  RVE  containing  diagonal  cracks  with  twice  the  crack  length  as  the  other  orientations 
(Figure  9e)  produced  slightly  more  accurate  results. 


Table  2.  Strain  Energy  Errors 


Technique 

Strain  Energy  (MPa) 

Error  (%) 

Fine  Mesh 

882.4 

- 

Figure  9a 

213.2 

75.8 

Figure  9b 

213.3 

75.8 

Figure  9c 

214.0 

75.7 

Figure  9d 

214.0 

75.7 

Figure  9e 

216.8 

75.4 

We  next  examine  both  the  homogenization  and  the  enrichment  techniques  using  900  quadratic  2nd  order 
elements,  again  examiningthestrainenergy  ofthemodel.Table3showsthatthehomogenization  method 
using  the  damaged  RVE  does  not  show  much  improvement  in  the  strain  energy  calculations,  with  lessthan 
a  %  difference  between  the  two  simulations.  Using  the  enrichment  method  without  any  damaged  RVEs 
still  reduced  the  error  by  6%  overthe  homogenization  approaches.  When  usingthe  enrichment  method 
along  with  the  damaged  RVE  reduced  the  error  to  60.4%,  a  15%  improvement  overthe  homogenization 
method. 


Table  3.  Strain  Energy  Errors 


Technique 

Strain  Energy  (MPa) 

Error  (%) 

Fine  Mesh 

882.4 

- 

Flomogenization  w/o 
Damaged  RVEs 

213.1 

75.8 

Flomogenization  with 
Damaged  RVEs 

216.8 

75.4 

Enrichment  w/o 
Damaged  RVEs 

273.5 

69.0 

Enrichment  with 
Damaged  RVEs 

349.7 

60.4 
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CONCLUSIONS 


In  this  paper  we  have  derived  the  enrichment  functions  for  microdamage  in  a  heterogeneous  medium. 
This  method  uses  a  structural  based  enrichment  technique,  allowing  macroscale  computations  to  be 
performed  with  the  microstructural  features  explicitly  considered.  In  section  three  we  have  demonstrated 
that  using  the  enriched  approach  with  damaged  RVEs  reduces  the  error,  for  the  problem  shown,  by  15%. 
We  are  continuing  this  research  for  dynamic  problems  and  are  expectingto  compare  the  results  with 
experimental  data. 


REFERENCES 

1.  Nemat-Nasser,  S  &  Hori,  M.  Micromechanics:  Overall  Properties  of  Heterogeneous  Materials,  North- 
Holland,  London,  1993. 

2.  Bakhalov,  N.  &  Panasenko,  G.  Homogenization:  Averaging  process  in  periodic  media.  Dordrecht:  Kluer: 
Academic  Publishers,  1989. 

3.  Fish,  J.  &  Yuan,  Z.  "Multiscale  enrichment  based  on  partition  of  unity."  International  Journal  for 
Numerical  Methods  in  Engineering  24  (2005):  1341-1359. 

4.  Macri,  M.  &  De,  S.  "An  octree  partition  of  unity  method  (OctPUM)  with  enrichments  for  multiscale 
modeling  of  heterogeneous  media."  Computers  &  Structures  86  (2008):  780-795. 

5.  Macri,  M.  &  Littlefield  A.  "Enrichment  Based  Multi-scale  Modeling  of  Composite  Materials  undergoing 
Thermo-Stress",  International  Journal  for  Numerical  Methods  in  Engineering  93  (2013):  1147-1169. 


12 


Approved  for  public  release;  distribution  is  unlimited. 


